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Abstract 

A method for generating finite-dimensional approximations to the solutions of optimal control prob- 
lems is introduced. By employing a description of the dynamical system in terms of its attainable sets 
in favor of using differential equations, the controls are completely eliminated from the system model. 
Besides reducing the dimensionality of the discretized problem compared to state-of-the-art collocation 
methods, this approach also alleviates the search for initial guesses from where standard gradient search 
methods are able to converge. The mechanics of the new method are illustrated on a simple double 
integrator problem. The performance of the new algorithm is demonstrated on a 1-D rocket ascent 
problem (” Goddard Problem”) in presence of a dynamic pressure constraint. 

Introduction 

The precise solution of an optimal control problem via Pontryagin’s Minimum Principle involves the 
numerical treatment of highly nonlinear multipoint boundary value problems (BVPs). The structure of 
these BVPs depends on the sequence in which the optimal control switches between singular/nonsingular 
and constrained/unconstrained arcs, and is not known to the analyst in advance. Additionally, these 
BVPs involve artificial costates that have little physical meaning, so that reasonable initial guesses for 
gradient search methods may be hard to fmd. 

For these reasons, rapid trajectory prototyping is usually attempted by applying direct optimization 
techniques to some type of discretized problem formulation. This approach leads to the numerical task 
of solving nonlinear programming problems. The performance of the optimization algorithm involved 
and the precision of the obtained solutions depends strongly on the chosen problem discretization and 
on the dimension of the associated parameter space. 

The currently most successful approaches are based on collocation methods [1], [2], as implemented 
in the OTIS program (Optimal Trajectories by Implicit Simulation) [3], [4]. The algorithm introduced 
in the present paper is very similar in its structure to the OTIS approach. However, the derivation 
is very different and involves concepts such as the hodograph space and differential inclusions. The 
advantage of the new approach lies in the fact that it is completely devoid of controls and, hence, 
requires a lower-dimensional parameter space than the OTIS approach. Furthermore, the absence of 
controls reduces the number of initial guesses required by nonlinear programming methods. 


* Senior Project Engineer, Analytical Mechanics Associates, Inc., 17 Research Drive, Hampton, VA 23666, working 
under contract at the Spacecraft Controls Branch, NASA-LARC, Member AIAA. 
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Problem Formulation 
in terms of Differential Equations 

We consider optimal control problems of the following general form: 



min $(i(0),ar(l)) 

(1) 

subject to the equations of motion 

i(t) = /(*(<), u (t)) V < G [0, 1], 

(2) 

the boundary conditions 

¥(*(0),*(1))=0, 

(3) 

and the control constraints 

g(x(t),u(t )) = 0 V t G [0, 1], 

0) 


h(x(t),u(t)) < 0 V< € [0,1]. 

(5) 


Here, t G R, x(t) 6 R n , and u(t) € R m are time, state vector and control vector, respectively. The 
functions $ : R 2 ” — R, / : R” +m -> R n , $ : R 2n - R% 5 < 2 n, and g : R n+m - R k ',h: K n+m -* 
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R*h are assumed to be sufficiently smooth w.r.t. their arguments of whatever order is required in this 
paper. (PWC[ 0, l]) m denotes the set of all piecewise continuous functions mapping the interval [0, 1] 
into R m . 


The Hodograph 

For fixed states x, the hodograph S(x) is defined as the set of all possible state rates x that can be 
achieved by varying the controls within their allowed bounds. Given the state equations (2) and the 
control constraints (4), (5), we can write 

S(: r) = {x e R" | x = f(x,u ), u e Jl(ar)} , (6) 

where f l(x) is the set of all admissible controls u 6 R m , i.e. 

Q(x) = {a£ R m | g(x,u) = 0, h(x,u) < 0}. (7) 

In the formulation (7), the controls u € R m can be regarded as a mere instrument for parameterizing 
the hodograph (6). In fact, the optimal state history and the optimal cost associated with problem (1) - 
(5) are not affected if the control vector u € R m and the control constraints (4), (5) are replaced by 
any other set of variables to parameterize the admissible state rates as long as the resulting hodograph 
remains unchanged. 

We now assume that there are smooth functions p : R 2 " — > R ,p , q ' R 2n — *• R-S such that the 
hodograph S(x) defined in (6), (7) can be rewritten as 

S(x) = {xe R n | p(M) = 0, q{x,x) < 0} . (8) 

To guarantee the existence of such functions p and q, we may replace S' ( x ) by its convex hull. This is 
sometimes called “relaxing the problem”. If the solution to the relaxed problem has its state rates always 
operating in the domain of the original nonconvex hodograph, then obviously the problem relaxation 
did not have any effect on the solution. In this case, the solutions to the relaxed and the unrelaxed 
problem are the same. If the solution to the relaxed problem has state rates operating outside the 
original nonconvex hodograph, then this indicates chattering control [5] and a solution to the unrelaxed 
problem does not exist. 

Hence, to avoid the problem of nonexistence of a solution to the original problem (1) - (5), we may, 
without loss of generality, assume that the hodograph defined in (6), (7) is convex. 

The aim of the sections below is to exploit (8) to introduce a problem formulation that is completely 
devoid of controls, employing only the information condensed in the hodograph. Using the concept 
of differential inclusions [6], [7], the evolution of the dynamical system can be described completely in 
terms of states and their sets of attainability. 

Sets of Attainability 

Given a starting time t 0 , an initial state x(t 0 ) = x 0 , and a final time <i, the set of attainability 
A'(f 0 ,zo;<i) is defined to consist of all points x 6 R" to which the state vector x(t) can be steered at 
time through an admissible control n(f-)t€[Wj] - Here, an admissible control is any function of time 
u(t) £ PWC\ 0, 1] restricted to the subinterval [/ 0 ,*i], that satisfies the control constraints (4), (5). The 
dependence of the attainable set on the right-hand side of the state equations /, and on the control 
constraints g , h, is suppressed in the nomenclature K(to,xo, fi). 

Now, let At be a small time step initiated at time to and let t\ = to + At. Then, a first-order 
approximation K(to,Xo]ti) to the attainable set Ii(to,xo\t\) can be given by 

K(t 0 ,x 0 ] <i) = G R n | x = x 0 + At ■ S(x 0 )} . (9) 

Here, we used the simplifying notation 

{At-S(x 0 )}:={&t-s\s£S{x 0 )}, (10) 
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Wi 


and S denotes the hodograph defined in (8). The notion of first, order approximation is understood in 
the sense that for fixed A t max , there is a real number M(At max ) > 0 such that for all 0 < A t < A t max 
each element of the attainable set K(t 0 ,x 0 ,ti) can be approximated to first order by some element in 
A'(/ 0 >£o?<i)« This means, for each element x 2 £ K(t 0 ,x 0 ,ti) there is an element x 2 £ K(t 0 ,x 0 ,ti) such 
that \\xi - i\ 1 1 2 < M • At 2 , 


A New Numerical Approach 

Let 0 = = 1 be a user-chosen subdivision of the interval [0, 1]. For simplicity, let 

the nodes be distributed equidistantly, i.e. 


A generalization to arbitrary subdivisions is straightforward. Then let x 0 ,xi, € R n be approxi- 
mations to the states at time respectively, and define 



X = [x 0 |xi| ■ -|*jv] € R (iV+1) ' n . 

(12) 

Problem (1) - (4) can now 

be rewritten in the form 



min $(£o? xaj) 

X € R(*+o-n v 7 

(13) 

subject to the constraints 

V(x 0 ,x N ) - 0 

(14) 

and 

*i+i E K(ti,Xi,U + i), i — 0, N — 1. 

(15) 

By substituting the attainable set x, - , i;+i) by its first-order approximation K(ti, x,-, <»+i) as 

by (9), (15) leads to 

x,-+i E {x € R" | x = Xi + At • 5(x,)} . 

given 

(16) 


Employing the assumption that the hodograph S(x) defined in (6) and (7) can be expressed in the 
form (8), it is clear that (16), after approximating by 


Xi = 


® t+1 ® i 


At 


can be substituted equivalently by the conditions 


= 0 

q(xi,Xi ) < 0 


* = 0,..,N- 1. 


(17) 


(18) 


In summary, it is proposed to obtain the approximate values of the optimal state vectors x % at 
times ti, i = by solving the nonlinear parameter optimization problem (13) subject to the 

constraints (14) and (18). The values of the state variables Xi at the node points Z,, i = 0, ,.,iV 
are the parameters that have to be found such that the performance index (13) is minimized. No 
controls are involved explicitly. The minimization is subject to the boundary conditions (14) and 
the interior constraints (18). For every state x* that is picked at time the set of attainability 
K(t„Xi\t) y i.e. the set of coordinates to which the state vector can be steered after time keeps 
ballooning as time progresses. The conditions (18) enforce that the state x; +2 lies within the (first-order 
approximation to) the attainable set A'(*,, x t ; / t+i ) (see figure 1). Hence, conditions (18) represent a 
numerical implementation of the differential inclusion concept. 
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Extension 1: 

A Higher Order Approximation 

A more precise discretization can be obtained if condition (16) is substituted by 

x, + i € {x £ R n I X - Xi + At • S(x,)} , (19) 


where 


Xi 


x, + x; + i 
2 


In complete analogy to (18) this leads to 


p(x„x t ) = 0 

q(ii,X{) < 0 


i = 0,..,JV- 1. 


A derivation of the order of this approximation is not given in the present paper. 


( 20 ) 


( 21 ) 


Extension 2: 

Nonequidistant Subdivision 

For practical applications it may be useful to place the nodes ti nonequidistantly rather than as 
defined in (11). For instance, if preliminary results obtained by equidistant node placement suggest 
rapid state transitions in some domain of the time interval, then it is advisable to rerun the problem 
with the same number of nodes, placed more densely in the areas of rapid state transitions and more 
scarcely in areas with sluggish state rates. Without increasing the number of parameters used to 
represent the discretized problem, this can significantly improve the precision of the result. 

Formally, nonequidistant node placement does not complicate the discretized problem formula- 
tion (13), (14), (16). 


Extension 3: 
State Constraints 

State equality or inequality constraints of the general form 

v(x) = 0, 
w(x) < 0, 


( 22 ) 


usually represent a significant complication of the optimal control problem (1) - (5). For the discretiza- 
tion proposed in this paper, constraints (22) are virtually trivial. By enforcing pointwise satisfaction 
of (22) we obtain the additional conditions 


v(xi) = 0 , 

u>(x,) < 0, 


i — 0, .., N. 


(23) 


Then, the suboptimal solution to problem (1) - (5), and (22) is obtained by simply adding con- 
straints (23) to the nonlinear programming problem of (13), (14), and (18). In contrast to optimal 
control approaches based on the Minimum Principle [8] , [5], [9], the user need not provide any guesses 
of the optimal switching structure. 


Extension 4: 

Analytical Derivatives 

The numerical approach proposed in the present paper does not require explicit integration of the 
equations of motion. Instead, a number of equality and inequality constraints is imposed on any pair 
of neighboring states. As a consequence, the partial derivatives of the cost gradients and the constraint 
gradients, required by any Newton type method to solve the Kuhn- Tucker conditions associated with 
problem (13), (14), (18), (23), can be calculated analytically as long as the functional dependencies of 
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p, p , q, v, w on their arguments are known. With this rather easy access to analytical partial 
derivatives of the cost and constraint functions associated with the discretized optimization problem, 
the expensive evaluation of partied derivatives through finite differences can be eliminated. It can also 
be expected that analytical differentiation provides higher precision, which may be a deciding factor in 
case of a badly conditioned problem. 


Extension 5: 

Explicit Time Dependence 

In the problem formulation (1) - (5) and (22), no explicit time dependence of the describing functions 
$,$,/, g , h, v, and w is assumed. This does not represent a serious restriction. Explicit dependence 
of the right-hand side of the state equations / on time t, for example, can be transformed away by 
introducing the additional state equation and initial condition 


f = 1, 
r(0) = 0, 


(24) 


thus providing a state carrying the value of the current time. Variable final time problems can be dealt 
with by introducing the additional state T through 

T = 0 (25) 

and multiplying the right-hand side / associated with all other states with T . 

These techniques are very common and well known, and they can be applied to transform general 
optimal control problems to problems of the form (1) - (5) and (22). However, for each additional 
state introduced, the number of parameters in the discretized formulation (13), (14), (18), and (23), 
increases by N + 1, where N is the user chosen number of discretization nodes. However, in this 
discretized formulation it is not necessary to explicitly carry along conditions of the type (24) and 
(25). Through analytical integration, conditions (24) can be eliminated completely, and, in case of 
condition (25), the unknown constant of integration gives rise to a single scalar parameter that has 
to be added to the optimization parameters (12). In general, to keep the number of parameters in 
the nonlinear programming problem (13), (14), (18), and (23) small, it is advisable to customize the 
numerical approach by using analytical integration whenever possible. The implications for the analysis 
outlined in the sections above are rather straightforward and the numerical benefits may be worth the 
extra effort . 


Example 1: 

Double Integrator Problem 

As an academic example to demonstrate the general procedures required by the new approach, we 


consider the following problem: 

' min — x(l) 

u€PtVC[0,l] 

(26) 

subject to the equations of motion 

X — V, 

(27) 


V — u, 

the initial and final conditions 

x(0) = 0 , 

u(0) = 0, v(l) = 0, 

(28) 

and the control constraints 

- 1 < u < 1. 

(29) 


The optimal control solution to this problem is of a bang-bang type. The associated state and control 
time histories are given in figures 2, 3. 
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To apply the proposed algorithm, we first choose an integer N and define N + 1 (equidistantly 
placed) nodes 

U = i = o A. (30) 

Then the values of the states at the nodes i = 0 ,..,]^, are obtained from solving the 

constrained parameter optimization problem 


subject to the constraints 


and 


Here, X{, u,-, x,, Vi, i = 0 , N 


min 
[(*<> w) i=0 


X 0 = 0, 

Vo = 0, 

x,- — Vi = 0 

Vi - 1 < 0 

—ii — 1 < 0 



vn = 0, 

1 i = 0,..,N- 1. 


1 are defined by 


Vi 


2 ’ 
TJi+l +V| 

2 ’ 




3?i4- 1 
At ’ 

Vj+l -Vi 
At > 


(31) 


(32) 


(33) 


(34) 


with 

At = (35) 

and give approximate values for states and state rates in between nodes. Conditions (32) represent the 
initial/final conditions (28) in the discretized form, and conditions (33) replace the description, (27) and 
(29), of the underlying dynamical system. For the derivation of conditions (33), the hodograph defined 
in (6) and (7) 

5(x, v) = |[x, v] £ R 2 1 x = v y v = u, uefl}, (36) 

fi(x, v) = {u G R | — 1 1}. (37) 

is rewritten in the general form (8) 


5(x, v) = {[ir, v] G R 2 
|i — v ~ 0, v — 1 < 0, — v — 1 < 0, } 


(38) 


Then the conditions i-u = 0, v— 1<0, -v - 1 < 0 in (38) are, loosely speaking, evaluated in between 
subsequent nodes to yield (33). 

It is well known that the optimal solution to problem (26) - (29) is of a bang-bang nature with 
u(t) = Tl for 0 < t < 0.5 and u(t) = —1 for 0.5 < t < 1 (see figures 2 and 3, or figures 4 and 5). From 
the linearity of the state equations (27), it follows that the discretized solution, i.e. the solution to (31) - 
(35), will match the optimal solution perfectly as long as the control associated with the optimal solution 
is constant throughout each discretization interval. Noting that the optimal solution has a switching 
point only at time t = 0.5 and is identically constant elsewhere, it is clear that the discretized solution 
is identical to the optimal solution if and only if N is an even integer (odd number of nodes). This 
observation is also confirmed by the numerical results (see figures 2 and 3 for odd numbers of nodes, 
and figures 4 and 5 for even numbers of nodes). All numerical results were obtained by employing the 
nonlinear programming code, NPSOL [10], to solve the nonlinear programming problem (31) - (33). 
The generation of decent initial guesses for the parameters (z;, Vi) i=0 N , required by NPSOL, was no 
problem. Convergence was always achieved in no more than three iterations even if the initial guesses 
were chosen many orders of magnitude off the respective optimal values. 
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Example 2: 

1-D Rocket Ascent (“Goddard Problem”) 

As a nontrivial problem to demonstrate the performance of the new algorithm, we consider the 
problem of maximizing the final altitude for a sounding rocket ascending vertically under the influence 
of atmospheric drag and an inverse-square gravitational field. The states are radial distance r, velocity 
v, and mass m. The thrust magnitude T is the only control and is subject to fixed bounds 0 < T < T max 
(control constraints) and a dynamic pressure limit q < q m ax (state constraint). 

In nondimensional form the problem is given as follows: 


min - r(tj) 

subject to the equations of motion 

r — v 

T - D 1 

V = y 

m r* 

T 

fn = 

c 


(39) 


(40) 


the control constraint 

T G [0, T max ] 

(41) 

the boundary conditions 

a) r(0) — 1 

d) r(tj) to be maximized 


b) v{0) = 0 

e) v(tj) free 

(42) 

Sis 

o 

II 

t— * 

/) m (tj) = mj 


and the state inequality constraint 

i 


v - 

/ ^ Qmax x n 

V-i £0 - 

(43) 


With dynamic pressure q and air density p given by 


9 = 



3 


respectively, it is clear that the “speed limit” (43) is equivalent to a dynamic pressure limit q-q ma x < 0. 
The aerodynamic drag D is given by 

D = q C D A. 


The constants Cd, A, p 0 , 0 denote drag coefficient, cross-sectional area, air density at ground level, 
and exponential decay rate of air density with altitude, respectively. The constants c, T max , mj used 
in (40), (41), (42) denote the exhaust velocity, the maximum available thrust, and the final mass of 
the vehicle (after all the fuel is burned), respectively. The nondimensional values used for numerical 
calculations are as follows: 

Cd 

Po ■ A 

P 

c 

Tmax 

m j 

A precise treatment of the problem above employing optimal control techniques is presented in [11]. 
For q max = oo, the time histories of the optimal states and the associated dynamic pressure are given 
in figures 6 - 9. For q max = 10, the results are shown in figures 10 - 13. 


0.05, 

12400, 

500, 

0.5, 

3.5, 

0 . 6 . 


(44) 
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To apply the numerical techniques introduced in the previous sections, we first choose an integer N 
and define N + 1 (equidistantly placed) nodes 


f'i — > i — o, N. 


(45) 


Then the values of the states at the nodes U, i = 0 are obtained from solving the 

constrained parameter optimization problem 

min — rjv (46) 


subject to the constraints 


7*0 = 

1, 


^0 = 

o, 


mo = 

= 1, 

m 

f% - Vi 


0 

tj. 4. zA i±R 

^ mi 

= 

0 

rhi 

< 

0 

-fhi - 

< 

0 


and 


v l - 


Poe 


■l gmaiE - <0 i = 0 N 

0(i -fi) - ’ 


Here, for all states x £ (r, v, m} , the quantities Xi, Ij, i = 0 , N - 1 are defined by 


_ £jii±£i. k..— 

— 2 » ~ At ’ 


(47) 


(48) 


(49) 


(50) 


with 

(51) 

and give approximate values for states and state rates in between nodes. Conditions (47) represent the 
initial/final conditions (42) in the discretized form, conditions (48) replace the description (40) and (41) 
of the underlying dynamical system, and conditions (49) enforce the state inequality constraint (43). 
Note that in contrast to the example problem 1, the final time tj is free here and appears as an additional 
parameter in (46) - (51). 

A first numerical solution is generated for the simple case N = 2. Using the rough initial guesses 

(ro n r 2 ) = (11 1) 

(v 0 Vl v 2 ) = (0 0 0) / 52 ) 

(mo mi m 2 ) = (1 0.8 0.6) 

h = 0.1 

the nonlinear programming problem (46) - (49) converges after less than 10 iterations with the code 
NPSOL [10]. Initial guesses for cases with N > 2 were generated by linearly interpolating the results 
obtained with N = 2. For q max = oo, figures 6 - 9 show the states r, v, m, and the dynamic pressure q 
versus time t, respectively. Figures 10 - 13 show an active state constraint case with q max = 10. 
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Summary and Conclusions 

A method for generating approximate solutions to optimal control problems has been introduced 
in this paper. By employing the concepts of attainable sets and differential inclusions, a numerical 
representation of the dynamical system is achieved that is completely devoid of controls. This leads to a 
discretized problem formulation of relatively low dimensionality. The absence of fast “moving” controls 
also improves the convergence properties and enhances robustness. 

The new method is illustrated on a detailed treatment of a simple double integrator problem. The 
performance of the algorithm is demonstrated on a 1-D rocket ascent problem (“Goddard Problem”) 
with and without an active dynamic pressure limit. 
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Figure 1: Schematic description of the differential inclusion approach 
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Figure 2: Example 1: State x vs time t for several odd node numbers 
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Figure 4: Example 1: State x vs time t for several even node numbers 
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Figure 5: Example 1: State v vs time t for several even node numbers 
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Figure 6: Example 2: State r vs time t for q max 
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Figure 10: Example 2: State r vs time t for q max 
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Figure 12: Example 2: State m vs time t for q max = 10 
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Figure 13: Example 2: Dynamic pressure q vs time t for q ma x = 10 
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